﻿//TODO: implement and Test, Wiener

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;
using System.Drawing.Imaging;

using LowLevelGraphics.Types;

namespace LowLevelGraphics.Filter
{   
    //TODO
    /// <summary>
    /// Restoration filter.
    /// Tries to recreate some lost informations by a certain frequency transform.
    /// </summary>
    public class Wiener : BaseImageFilter
    {
        /// <summary>
        /// Initializes a new instance of the <see cref="Wiener"/> class.
        /// </summary>
        public Wiener()
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="Wiener"/> class.
        /// </summary>
        /// <param name="_wiener">The _wiener.</param>
        internal Wiener(Wiener _wiener)
        {
        }

        /// <summary>
        /// Widthes the bytes.
        /// </summary>
        /// <param name="_nBits">The _n bits.</param>
        /// <returns></returns>
        private int WIDTHBYTES(int _nBits)
        {
            return ((((_nBits) + 31) >> 5) << 2);
        }

        /// <summary>
        /// Execute this Filter
        /// A Debug Statement is written if a concrete Filter is called
        /// </summary>
        /// <param name="_bitmap"></param>
        /// <returns></returns>
        public override UnsafeBitmap Execute(UnsafeBitmap _bitmap)
        {
            UnsafeBitmap bitmap = _bitmap;

            int lHeight = bitmap.Height;
            int lWidth = bitmap.Width;
            //
            string lpSrc;
            //
            long i;
            long j;
            //
            byte pixel;
            //
            long lLineBytes;
            //
            double[] fftSrc, fftKernel, fftNoise;
            double a, b, c, d, e, f, multi;

            //
            ulong[] nn = new ulong[3];
            //
            double MaxNum;

            //
            lLineBytes = WIDTHBYTES(lWidth * 8);

            ////if we got not the correct size here, the pic should be resized then.
            double dPower = Math.Log((double)lLineBytes) / Math.Log(2.0);
            //if (dPower != (int)dPower)
            //{
            //    return _bitmap;
            //}

            //dPower = Math.Log((double)lHeight) / Math.Log(2.0);
            //if (dPower != (int)dPower)
            //{
            //    return _bitmap;
            //}

            fftSrc = new double[lHeight * lLineBytes * 2 + 1];
            fftKernel = new double[lHeight * lLineBytes * 2 + 1];
            fftNoise = new double[lHeight * lLineBytes * 2 + 1];

            nn[1] = (ulong)lHeight;
            nn[2] = (ulong)lLineBytes;
            
            for (j = 0; j < lHeight; j++)
            {
                for (i = 0; i < lLineBytes; i++)
                {
                    
                    lpSrc = string.Empty; //(char*)lpDIBBits + lLineBytes * j + i;

                    pixel = 0;//(uchar) *lpSrc;

                    fftSrc[(2 * lLineBytes) * j + 2 * i + 1] = (double)pixel;
                    fftSrc[(2 * lLineBytes) * j + 2 * i + 2] = 0.0;

                    if (i < 5 && j == 0)
                    {
                        fftKernel[(2 * lLineBytes) * j + 2 * i + 1] = 1 / 5.0;
                    }
                    else
                    {
                        fftKernel[(2 * lLineBytes) * j + 2 * i + 1] = 0.0;
                    }
                    fftKernel[(2 * lLineBytes) * j + 2 * i + 2] = 0.0;
                    if (i + j == ((int)((i + j) / 8)) * 8)
                    {
                        fftNoise[(2 * lLineBytes) * j + 2 * i + 1] = -16.0;
                    }
                    else
                    {
                        fftNoise[(2 * lLineBytes) * j + 2 * i + 1] = 0.0;
                    }
                    fftNoise[(2 * lLineBytes) * j + 2 * i + 2] = 0.0;
                }
            }

            //srand((unsigned)time(NULL));
            //
            Fourn(fftSrc, nn, 2, 1);
            ////
            Fourn(fftKernel, nn, 2, 1);
            ////
            Fourn(fftNoise, nn, 2, 1);

            for (i = 1; i < lHeight * lLineBytes * 2; i += 2)
            {
                a = fftSrc[i];
                b = fftSrc[i + 1];
                c = fftKernel[i];
                d = fftKernel[i + 1];
                e = fftNoise[i];
                f = fftNoise[i + 1];
                multi = (a * a + b * b) / (a * a + b * b - e * e - f * f);
                if (c * c + d * d > 1e-3)
                {
                    fftSrc[i] = (a * c + b * d) / (c * c + d * d) / multi;
                    fftSrc[i + 1] = (b * c - a * d) / (c * c + d * d) / multi;
                }
            }

            //
            //fourn(fftSrc, nn, 2, -1);

            //
            MaxNum = 0.0;
            for (j = 0; j < lHeight; j++)
            {
                for (i = 0; i < lLineBytes; i++)
                {
                    fftSrc[(2 * lLineBytes) * j + 2 * i + 1] =
                        Math.Sqrt(fftSrc[(2 * lLineBytes) * j + 2 * i + 1] * fftSrc[(2 * lLineBytes) * j + 2 * i + 1] + fftSrc[(2 * lLineBytes) * j + 2 * i + 2] * fftSrc[(2 * lLineBytes) * j + 2 * i + 2]);
                    if (MaxNum < fftSrc[(2 * lLineBytes) * j + 2 * i + 1])
                        MaxNum = fftSrc[(2 * lLineBytes) * j + 2 * i + 1];
                }
            }

            //
            for (j = 0; j < lHeight; j++)
            {
                for (i = 0; i < lLineBytes; i++)
                {
                    // 	
                   // lpSrc = (char*)lpDIBBits + lLineBytes * j + i;

                    //*lpSrc = (uchar)(fftSrc[(2 * lLineBytes) * j + 2 * i + 1] * 255.0 / MaxNum);
                }
            }
            //
            return _bitmap;
        }

        /// <summary>
        /// Fourns the specified data.
        /// </summary>
        /// <param name="data">The data.</param>
        /// <param name="nn">The nn.</param>
        /// <param name="ndim">The ndim.</param>
        /// <param name="isign">The isign.</param>
        /// <returns></returns>
        protected bool Fourn(double[] data, ulong[] nn, int ndim, int isign)
        {
            int idim;
            ulong i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
            ulong ibit, k1, k2, n, nprev, nrem, ntot;
            double tempi, tempr;
            double theta, wi, wpi, wpr, wr, wtemp;

            for (ntot = 1, idim = 1; idim <= ndim; idim++)
                ntot *= nn[idim];
            nprev = 1;
            for (idim = ndim; idim >= 1; idim--)
            {
                n = nn[idim];
                nrem = ntot / (n * nprev);
                ip1 = nprev << 1;
                ip2 = ip1 * n;
                ip3 = ip2 * nrem;
                i2rev = 1;
                for (i2 = 1; i2 <= ip2; i2 += ip1)
                {
                    if (i2 < i2rev)
                    {
                        for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2)
                        {
                            for (i3 = i1; i3 <= ip3; i3 += ip2)
                            {
                                i3rev = i2rev + i3 - i2;
                                Sorter.Swap<double>(ref data[i3], ref data[i3rev]);
                                Sorter.Swap<double>(ref data[i3 + 1], ref data[i3rev + 1]);
                            }
                        }
                    }
                    ibit = ip2 >> 1;
                    while (ibit >= ip1 && i2rev > ibit)
                    {
                        i2rev -= ibit;
                        ibit >>= 1;
                    }
                    i2rev += ibit;
                }
                ifp1 = ip1;
                while (ifp1 < ip2)
                {
                    ifp2 = ifp1 << 1;
                    theta = isign * 6.28318530717959 / (ifp2 / ip1);
                    wtemp = Math.Sin(0.5 * theta);
                    wpr = -2.0 * wtemp * wtemp;
                    wpi = Math.Sin(theta);
                    wr = 1.0;
                    wi = 0.0;
                    for (i3 = 1; i3 <= ifp1; i3 += ip1)
                    {
                        for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2)
                        {
                            for (i2 = i1; i2 <= ip3; i2 += ifp2)
                            {
                                k1 = i2;
                                k2 = k1 + ifp1;
                                tempr = wr * data[k2] - wi * data[k2 + 1];
                                tempi = wr * data[k2 + 1] + wi * data[k2];
                                data[k2] = data[k1] - tempr;
                                data[k2 + 1] = data[k1 + 1] - tempi;
                                data[k1] += tempr;
                                data[k1 + 1] += tempi;
                            }
                        }
                        wr = (wtemp = wr) * wpr - wi * wpi + wr;
                        wi = wi * wpr + wtemp * wpi + wi;
                    }
                    ifp1 = ifp2;
                }
                nprev *= n;
            }
            return true;
        }

        /// <summary>
        /// This is used to clone a filter, each filter must implement it's own clone functionality
        /// </summary>
        /// <returns></returns>
        public override object Clone()
        {
            return new Wiener(this);
        }
    }
}
